import numpy as np
from scipy.stats import norm


class BSMFDM:
    """Black Scholes Model
    Example:

    K = 100;  r = 0.1 ;   T = 1 ;   sigma = 0.3 ;   S = np.arange(60,140,0.1)

    https://www.codearmo.com/python-tutorial/options-trading-greeks-black-scholes

    Attributes:
        S : current asset price
        K: strike price of the option
        r: risk free rate
        T : time until option expiration
        σ: annualized volatility of the asset's returns
        q: continuously compounded dividend ,0 to ignore dividend

    Returns:
        Call: Call=S0*N(d1)*(e^-qT)-N(d2)K*(e^-rT)
        Put: Put=N(-d2)*K*(e^-rT)-N(-d1)*S0*(e^-qT)
        d1: d1=[ln(S/K)+(r-q+omega^2/2)*T]/[omega*(T^1/2)]
        d2: d2=d1-omega*(T^1/2)
    Raises:
        cpt err: every time new BSM instance to avoid grade cpt pollute attribution value.


    """

    def __init__(self, S, K, T, r, opt_ty="call", sigma=0.5, q=0.0):
        self.S = S
        self.K = K
        self.T = T
        self.r = r
        self.opt_ty = opt_ty
        self.sigma = sigma
        self.q = q
        self.ds = 1e-4

    @staticmethod
    def N(x):
        return norm.cdf(x)

    def d1(self):
        return (np.log(self.S / self.K) + (self.r - self.q + self.sigma ** 2 / 2) * self.T) \
               / (self.sigma * np.sqrt(self.T))

    def d2(self):
        return self.d1() - self.sigma * np.sqrt(self.T)

    def call_value(self):
        return self.S * np.exp(-self.q * self.T) * self.N(self.d1()) - \
               self.K * np.exp(-self.r * self.T) * self.N(self.d2())

    def put_value(self):
        return self.K * np.exp(-self.r * self.T) * self.N(-self.d2()) - \
               self.S * np.exp(-self.q * self.T) * self.N(-self.d1())

    def price(self):
        if self.opt_ty == 'call':
            return self.call_value()
        if self.opt_ty == 'put':
            return self.put_value()
        else:
            raise ValueError('Unrecognized type')

    def implied_sigma(self, market_price, one_cent=0.01, max_try=1000):
        """隐含波动率  牛顿逼近 本质就是梯度下降,从self.sigma起步搜索 默认0.5
            Args:
                option_ty: C for Call ; P for Put
                market_price: market price
                sigma0: start point for sigma
                max_try: how many loops to try
        """
        for i in range(max_try):
            _bs_price = self.price()
            diff = market_price - _bs_price
            if abs(diff) < one_cent:
                return self.sigma
            self.sigma += diff / self.vega()
        return self.sigma

    def delta(self):
        """基础价格对期权价格求导
           nd1 = self.N(self.d1())
        if option_ty == "c":
            return nd1 * np.exp(-self.q * self.T)
        else:
            return [nd1 - 1] * np.exp(-self.q * self.T)
        """
        p0 = self.price()
        self.S = self.S + self.ds
        p1 = self.price()
        return (p1 - p0) / self.ds

    def gamma(self):
        """ 基础价格对期权价格二阶导 self.ds过小可能不准
        """
        p0 = self.price()
        self.S = self.S + self.ds
        p1 = self.price()
        self.S = self.S + self.ds
        p2 = self.price()
        return (p2 - 2 * p1 + p0) / (self.ds * self.ds)

    def vega(self):
        """波动率对期权价格
        """
        p0 = self.price()
        self.sigma = self.sigma + self.ds
        p1 = self.price()
        return (p1 - p0) / self.ds

    def theta(self):
        """期限对期权价格
        到期时间不是时间减法,所以要加-
        """
        p0 = self.price()
        self.T = self.T + self.ds
        p1 = self.price()
        return -(p1 - p0) / self.ds

    def rho(self):
        """期限对利率价格
        """
        p0 = self.price()
        self.r = self.r + self.ds
        p1 = self.price()
        return (p1 - p0) / self.ds


if __name__ == '__main__':
    # calibrated by https://quantpie.co.uk/oup/oup_bsm_price_greeks.php
    S = 120
    K = 100
    r = 0.03
    opt_ty = "put"
    T = 0.5

    sigma = 0.5
    q = 0.1
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("price==", b.price())
    mprice = 9
    print("imply2", b.implied_sigma(mprice))
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("price==", b.price())
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("delta", b.delta())
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)

    print("gamma", b.gamma())
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("vega", b.vega())
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("theta", b.theta())
    b = BSMFDM(S, K, T, r, opt_ty, sigma, q)
    print("rho", b.rho())
